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Abstract 

The matched interface and boundary (MIB) method has a proven ability for delivering the second 
order accuracy in handling elliptic interface problems with arbitrarily complex interface geometries. 
However, its collocation formulation requires relatively high solution regularity. Finite volume method 
(FVM) has its merit in dealing with conservation law problems and its integral formulation works well 
with relatively low solution regularity. We propose an MIB-FVM to take the advantages of both MIB 
and FVM for solving elliptic interface problems. We construct the proposed method on Cartesian 
meshes with vertex-centered control volumes. A large number of numerical experiments are designed 
to validate the present method in both two dimensional (2D) and three dimensional (3D) domains. 

It is found that the proposed MIB-FVM achieves the second order convergence for elliptic interface 
problems with complex interface geometries in both Loo and L 2 norms. 
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1 Introduction 


Elliptic partial differential equations (PDEs) with discontinuous coefficients and singular source terms, 
commonly referred to as elliptic interface problems, occur in many applications, including fluid dynamics 
dn [2511311 in], material science [ 22 ] , electromagnetics [HIEIIEHIIM], biological systems, [HI EH H EH m 
and heat or mass transfer m- Since the pioneer work of Peskin in 1977 [45], lots of attention has been 
paid to this field in the past few decades 

An elliptic interface problem is formulated as elliptic PDEs defined on piecewise-smooth subdomains, 
which are coupled together via interface conditions, such as given jumps in solution and flux across the 
domain interface. Without the interface, a large variety of convergent numerical schemes are available 
for solving the elliptic PDEs, such as, finite difference method (EDM), finite element method (EEM), 
finite volume method (EVM), wavelet, radial basis functions, meshless and spectral methods. However, for 
elliptic interface problems, the direct application of the aforementioned schemes cannot yield a convergent 
solution. 

Peskin proposed the immersed boundary method (IBM) [T7l|30[|45] in order to simulate flow pattern of 
blood in the heart, in which he approximated the singular sources on the interface. Since his work, many 
numerical methods have been proposed and studied. In 1984, Mayo introduced a second order integral 
equation approach for Poisson’s equation and biharmonic equation on irregular domains [40l|4T], in which 
the solution is extended to a rectangular region by using Eredholm integral equations. A fast Poisson 
solver is utilized to solve the resulting Eredholm integral equations on a rectangular region. Continuous 
derivatives have been assumed to evaluate the discrete Laplacian. This method can deal with jump 
conditions of [u] 7 ^ 0 and [un] = 0 when the Green’s function is available. In 1994, a level set method in 
combination with the immersed boundary method [48] is proposed by Osher and his coworkers in order to 
compute solutions to incompressible two-phase flow. This level set method is easy to implement, although 
it is of low convergent rate. In the same year, immersed interface method (IIM) was proposed by Levique 
and Li lasiisiiii], in which the interface conditions are incorporated into the finite difference scheme 
near the interface to achieve second order accuracy based on a Taylor expansion in a local coordinate 
system. It is a second order interface scheme that does not smear the interface jumps and is used to solve 
elliptic and parabolic interface problems alike. The resulting linear system is sparse, but not symmetric 
or positive definite. Also, second order derivatives of the interface jumps are needed. The IIM is widely 
used in practice and its extensions can be seen in Refs. m^- Eor instance, in 2001, Li and Ito [34] 
constructed an IIM scheme with resulting linear matrix being diagonally dominate and its symmetric 
part being negative definite. 

In addition to the IIM, a large class of finite difference based numerical methods have also been 
proposed for solving elliptic interface problem. The crucial idea of these methods is to incorporate jump 
conditions into the difference stencils near the interface to maintain high order local truncation error 
using Taylor expansion. Numerical schemes based on finite element or finite volume are also developed 
[HElllSo]. Einite element schemes are usually constructed by modifying the finite element basis near 
the interface. A few methods on unfitted mesh have been introduced based on the discontinuous Galerkin 
method using interior penalty technique to handle jump and flux conditions 13 ng. The weak Galerkin 
finite element method has also been developed for elliptic interface problem [43]. Other interesting 
methods developed in recent decades include the ghost fluid method (GEM) by Eedkiw, Osher and 
coworkers [T2j [38] . The second order convergence of finite difference formalism of interface methods is 
proved in 2006 by Beale and Layton for smooth interface [4] , whereas convergence analysis of most elliptic 
interface schemes is yet to be done. 

In the past decade, we have dedicated ourselves to designing accurate and robust numerical schemes 
for solving elliptic interface problems. In 2006, the matched interface and boundary (MIB) method 
[63l|6ll[70l|69] was proposed, motivated by many practical needs, such as optical molecular imaging 
ig , nano-electronic devices, |g , vibration analysis of plates (Sg , wave propagation (HSl (Ml , geo dynamics 
[ 68 ] and electrostatic potential in proteins 113110111117]. One important feature of the MIB method 
is its extension of computational domains by using the so called fictitious values, an strategy developed 
in our earlier methods for handling boundaries ligisii and interfaces [ 66 ]. As such, standard central 
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finite difference schemes can be employed to discretize differential operators as if there were no interface. 
Another unique feature of the MIB method is to repeatedly enforce only the lowest order jump conditions 
to achieve higher order convergence, which is of critical importance for the robustness of the method to 
deal with arbitrarily complex interface geometries. Higher-order jump conditions must involve higher 
order derivatives and/or cross derivatives. Therefore, to approximate higher order derivatives and/or 
cross derivatives, one must utilize larger stencils, which is unstable for constructing high-order interface 
schemes and cumbersome for complex interface geometries. The other distinct feature of the MIB method 
is its dimensional splitting. To enforce a 2D or 3D interface jump condition, we divide the problem into 
multi ID ones and resolve a ID problem at a time, if possible. This approach enables us to come up 
with a systematic procedure to resolve high dimensional interface problems. Finally, based on high order 
Lagrange polynomials, the MIB method is of arbitrarily high order in principle. For example, MIB 
schemes up to 16th order accurate have been constructed for simple interface geometries ID and 2D 
domains [66] |70] , and sixth-order accurate MIB schemes have been developed for complex interfaces in 
2D [70] and 3D domains [61] |70] . Recently we have constructed an adaptively deformed mesh based MIB 
m and a Galerkin formulation of MIB [59] [56] to improve MIB’s capability of solving realistic problems. 
A comparison of the GFM, IIM and MIB methods can be found in in Refs. [701 [69]. 

The MIB method has been used in our earlier works for solving many scientific and engineering 
problems, such as Poisson-Boltzmann equation (PBE) [67] [60[ |T6l [7] for describing the electrostatic 
potential in proteins. To our best knowledge, the MIB method is the only method that has demonstrated 
the second order accuracy for solving Poisson or PB equation with realistic protein surfaces with geometric 
singularities Enmaiii and for solving Poisson equation with multiple material interfaces [58]. It can 
also be used to solve the Helmholtz equation for wave scattering and propagation in inhomogeneous 
media. A fourth order MIB scheme for the Helmholtz equation with arbitrarily curved interfaces has 
been proposed by Zhao [65] [64]. Another example is the geodynamics where the Navier-Stokes equations 
with discontinuous viscosity and density is to be solved. Zhou et al have developed a second order 
accurate MIB method to solve this problem on non-staggered Cartesian grids [68]. Furthermore, the 
elastic interface problem in both 2D and 3D domains with arbitrarily complex interface geometry was 
also addressed by the MIB method [51] [52]. In the past few years, the MIB method has also been applied 
to the optical molecular imaging [9], nano-electronic devices, [8] and vibration analysis of plates [62] . 

Due to the importance and complexity of interface problems, it is still urgent to develop methods 
that are more accurate, robust and require less computational cost. Finite volume method (FVM) is 
known for its ability to better conserve mass and fiux in conservation law problems. This property is 
fundamental for the simulation of many physical models, e.g., oil recovery simulation and computational 
fluid dynamics in general. Additionally, compared with finite difference method which takes a collocation 
formulation, FVM is able to deal with solutions with relatively lower regularity. In the present work, we 
propose the finite volume formulation of the MIB method (MIB-FVM). The motivation for the proposed 
MIB-FVM is to inherit the merits of both methods, namely, the capability of the MIB method in handling 
complex interface geometries and the advantage of FVM for conservation laws and low solution regularity. 

The rest of this paper is organized as the follows. In Section [^ the general theory of the MIB based 
finite volume formulation is briefly discussed. The theoretical formulation and the computational algo¬ 
rithms are given as well. In Section [^ we present the strategy for treating complex interface geometries. 
In Section [^ the proposed MIB-FVM is validated by benchmark tests, such as 6-petal flower and jigsaw 
puzzle-like shape in two dimensional (2D) space, and sphere, ellipsoid, cylinder, flower-based cylinder 
and torus in three dimensional (3D) space. Solutions of less regularity (i7^ continuous) are also tested 
for the spherical and ellipsoidal interface in 3D. This paper ends with a conclusion. 
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Figure 1: Illustration of a 3D interface problem at a cross section. The whole domain Q is divided by 
interface F into two subdomains, and . The mesh is represented by solid lines. The grid points 
are represented by blue dots. 

2 Matched interface and boundary - finite volume method (MIB- 
FVM) 

Let us consider an open bounded domain C with a given interface F, which separates the domain 
into two subdomains, and as illustrated in Figure [l] 

The boundary dQ and interfaces F may be non-smooth. The interface can be characterized by a 
piecewise smooth level-set function ip e such that F = {x|(^(x) = 0, Vx G As such, two subdomains 
can be given by 0+ = {x|(^(x) < 0,Vx G 0} and 0“ = {x|(^(x) > 0,Vx G The elliptic interface 
problem can be formulated as 

- V •/3(x)Vi^(x) = ^(x), Vx G (1) 

u = ^ 5 , Vx on dft, (2) 

where g{x.) is a piecewise continuous function, gij is the boundary value, and /3{x.) is a variable coefficient 
that is discontinuous across the interface F. As a result, two jump conditions are required to make the 
problem well posed 


[u] = — u = Vx on F (3) 

and 

[I3un] = I3+U+ - p-u- = Vx on T, (4) 

where and denote their limiting value from the side of the interface F, and u~^u~ and 

f3~ denote their limiting value from the side of the interface F. The derivatives and u~ are 
evaluated along the outer normal direction on the interface. Here ^(x) and T(x) is at least continuous. 
Equations §-@ define the elliptic interface problem to be solved in the present work. 
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Figure 2: Left chart: Illustration of the 3D control volume at a cross section that is perpendicular to the 
x-axis and intersect with the x-axis at the ith grid point. The yellow region denotes the control volume 
associated with the grid point (i, j, k). The grid points are denoted by blue dots. The mesh is represented 
by solid lines and the dual mesh is represented by the dashed lines. Right chart: Illustration of the 3D 
irregular domain at a cross section. The irregular domain is denoted by the yellow region. The grid 
points are denoted by blue dots. The mesh is represented by solid lines and the dual mesh is represented 
by dashed lines. 

2.1 Introduction to the MIB finite volume method 

To avoid the time-consuming mesh generation procedure, like the original MIB method, we employ 
the Cartesian mesh. We also employ the vertex-centered finite volume method, which means that we 
associates control volumes and unknowns to vertices. The control volumes are cubes whose centers are 
located at the grid points and whose sides intercept at the midpoints between grid points, which is 
illustrated by the left chart of Figure 

Like the original MIB method, the grid points and and control volumes are classified into regular 
and irregular types. A grid point is said to be irregular if the standard central finite difference (CFD) 
scheme at the grid point near the interface requires grid point (s) from the other side of the interface. 
In Figure [l] the interface F and the different domains ^1, (2+ and are depicted. As is shown in the 
figure, the complex interface inevitably cut through certain control volumes. A control volume is defined 
to be irregular if it is cut through by the interface, or stated differently, the vertices of the control volume 
locate on both side of the interface. The irregular domain E is defined to be the union of all irregular 
control volumes, as is shown in the right chart of Figure 

In the present MIB-FVM, we further define each control volume to be a positive or negative control 
volume according to where its center is, i.e., a control volume is said a positive one if its center grid 
point locates in the positive region similarly, a control volume is said a negative one if its center 
grid point locates in the negative region We also define and to be the union of all positive 
and negative control volumes, respectively, as is shown in Figure]^ Domain 0+ and can be viewed 
as the generalized 0+ and . First, domains and are disjoint by such definition; secondly, it 
can be easily seen that there is no inclusion relation between the origin subdomains or and the 
generalized subdomains or . 
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Figure 3: Illustration of the domain re-division at a cross section. Left chart: domain (in green), 
which is the union of all the positive control volumes; Right chart: domain (in yellow), which is the 
union of all the negative control volumes. 


To properly define the elliptic interface problem, we also need to generalize the discontinuous coefficient 
/3(x) and the source term ^(x) over the generalized domains. To this end, we define 


and 




Vx c f2+ n f2+ 


Vx c f2+ n f2+ 


/3 (x), 

Vx c nfi 


Vx c n 


(5) 


(6) 


where /d^(x) and /d^(x) are the smooth extensions of the coefficients /d^(x) and /3~(x) over domains 
and respectively. The superscript c means complement. 

Similarly, we define 


and 



5+(x), 

Vx c 0+ n n+ 


Vx c 0+ n n+ 


9 (x), 

Vx C 

9 ex (x) ; 

Vx c 0“ n 


(7) 

( 8 ) 


where and g~^ are the smooth extension of the functions g~^(x) and g~{x) over domains and , 
respectively. The superscript c denotes complement. 

By doing the generalization above, we can define the original interface problems on irregular control 
volumes. As the solution and its gradient can be discontinuous across the interface, to solve the interface 
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Figure 4: Illustration of the faces of control volume Vi. Left chart: Si and S 2 are two faces along the 
x-direction; middle chart: S 3 and S 4 are two faces along the ^-direction; right chart: S 5 and Sq are two 
faces along the ^-direction, rrii is the center of Si, where i = 1 , 2 ,..., 6 . 


problem, one needs to rigorously enforce the interface jump conditions. These conditions, interpreted by 
a set of equations, are utilized to determine fictitious values on generalized domains numerically. The 
detailed procedure of computing fictitious values is described in Section To achieve high order finite 
difference type of MIB schemes, the set of interface conditions are iteratively implemented to determine 
as many fictitious values as needed. This approach can also be used to create a high-order MIB-FVM. 


2.2 MIB finite volume formulation 

The MIB-FVM is based on the integral form rather than the differential form of Eq. 0. The 
integral form of the equation can be obtained by taking integral on both sides of the equation over any 
control volume Vi and using the divergence theorem on the left hand side 



(9) 


where dVi denotes the surface of the control volume V, and n denotes the unit outer normal vector of 
The above equation can be written as 



Puy,Pu^) ■ (ni,n2, 



g{x)dx, 


( 10 ) 


where (^ 1 , 712 ,^ 3 ) is the coordinate of the unit outer normal n. 

As is illustrated in Figure since the control volume Vi is a cube, dVi consists of 6 faces, denoted by 
dSi where i = 1, 2, ... 6 . For instance, dSi denotes the left face, whose unit outer normal is ( — 1,0, 0); dS 2 
denotes the right face, whose unit outer normal is ( 1 , 0 , 0 ), and so on. 

As a result, Eq. (10) can be further written as 






The above surface and volume integrals are further approximated by Gauss quadrature rules. In 
practice, we use the second order midpoint rule, thus the surface integrals on the left hand side are 
approximated by the function values at the face center, times the face area; the volume integral on the 
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Figure 5: Exemplary configurations for a regular control volume on the left and an irregular one on the 
right. Here i = 0,1, 2,6 in blue color denote the grid points, with Pq being the center of the control 
volume. Here rrii^i = 1,2, ...,6 in red color denote the centers of the corresponding faces. For the right 
chart, the green area represents the 2D interface in 3D that divides the control volume into two parts. 


right hand side is approximated by the function value at the center, times the volume of the control 
volume. Based on Figure and the midpoint rule, Eq. ( [IT] ) can be further written as 

-{Pux)rmdydz + {I3ux)m2dydz - {I3uy)msdxdz + {Puy)m^dxdz - {Puz)m^dxdy + {f3uz)medxdy = gidxdydz, 

( 12 ) 


where {^Ut)miP = x,y,z;i = 1,2, ...,6 denote the function value of {^Uf) at the ith midpoint rrii] Qi 
denote the function value of 5 f(x) at the center of the ith control volume Vf, and dx^dy and dz denote 
the size of partition in x, ^ and 2 : direction, respectively. 

For a regular control volume illustrated by the left chart of Figure the above partial derivatives at 
the midpoints can be evaluated by a 3-point interpolation scheme, namely 




• {wl,wl,wl) 

—5 ^i,j. 

\T 

,/c 5 '^i-\-l,j,k) 


« (/?)m2 

■ {wl,wl,wl) 


(/3Uy)m3 ^ 

^ {l3)m3 

' {wl,wl,wl) ■ 

—1,/c5 

\T 

k 5 ) 

{ldUy)m4 ^ 

■ 

■ {wl,wl,wl) ■ 


{PUz)m5 ^ 


■ {wl,wl,wl) 

,k — l’! ^i,j. 

\T 

k") ) 

{PUz)me ^ 

- iP)me 

■ iwl,wl,wl) 

{'^i,j,k — l’)'^i,j,k’)'^i,j,k-\-l) 


where wl^i = 1,2...,6; j = 1,2,3 are the interpolation weights, and denotes the function value at 
grid point {i^j^ k). 

For an irregular control volume illustrated by the right chart of Figure to maintain the same 
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convergence rate and accuracy, the above interpolations must be modified as 


(/^'^ai)mi ^ 



' {Ui- 


1 

\T 

) '^i+l,jf,/e ) 

{/3Ux)m2 ^ 

^ {^o)m2 

• ■ 

' {Ui- 


> 

\T 

> l^i+l,j,/e j 

(/3Uy)m3 ^ 

^ {Po)m3 

■ {wl,wl,wl) ■ 

(UiJ 


^id^k 5 

\T 

l^i,j)+l,/e ) 

(/3Uy)m4 ^ 

^ (/^o)m4 

■ {wl,wl,wl) • 


— 1,/e 5 

^id^k 5 

\T 

l^i,ji‘+l,/e j 

(/3Uz)m5 ^ 


• {wl,wl,wl) ■ 

i^i,j 

,/c-l; 

) '^id,k’ 

\T 

) l^i,j),/e+l ) 

(/3Uz)m6 ^ 

^ {^o)me 

■ {wl,wl,wl) ■ 


,/e —1; 


> ^id,k-\-l) 


where /3o denotes the generalization of (3 defined in Eqs. (§ and (|^, and u is defined as 


^x,y,z 


'^x,y,z if (^5 Hi z) is on the same side of the interface with {i^j, k) 
fx,y,z if is on different side of the interface with (i, j, k). 


( 14 ) 


Here fx,y,z represents the fictitious value at grid point (x,^, z)^ which is computed by MIB method. The 
detailed procedure of determining these fictitious values are described in the next section. 


3 Algorithms for determining fictitious values 


3.1 Simplification of interface jump conditions 

To solve the above linear equation system, we need to determine the involved fictitious values first. This 
issue is discussed in the present section. We first describe the interface jump conditions, which are needed 
at each intersecting point of control volume edges and the interface. 

As the interface normal direction varies along the interface, which is very troublesome from the 
computational point of view. We introduce a local coordinate system at each intersection point of 
the Cartesian mesh and the interface to treat different interface geometries systematically. At a specific 
intersection point, the local coordinate system is denoted by (^, C), where ^ denotes the normal direction 

and T] is on the x—y plane. The relation of the local coordinates and the Cartesian coordinates is described 
in the following 



where P is the transformation matrix 


P = 


sin (p cos 0 
— sin^ 

— cos (j) cos 0 


sin (j) sin 0 
cosO 

— cos (j) sin 0 



(15) 


(16) 


where 0 and (j) are the azimuth and zenith angles with respect to the normal direction n, respectively. 

Without generating higher order jump conditions, usually we differentiate Eq. (§ in two tangential 
directions in order to obtain two additional jump conditions [u^j] = — u~ ^ and [u(^] = u'^ — u'^. Thus, 

in the new coordinate system, the jump conditions in two tangential directions can be discretized as 


[ur^] = (— sin Ou'^ + cos OUy ) — (— sin Ou^ + cos OUy ) (17) 

= (— cos 0 cos ^14 J — co^(p^mOUy + sin^i^^) — (— co^<pco^0u~ — cos0sin^i4“ + ^m(pu~). (18) 

Similarly, in the new coordinate system, the jump condition Eq. @ can be discretized as 
[Pu^] = j3^ (sin 0 cos Ou^ + sin 0 sin OUy + cos (pu^) — /3~ (sin (p cos 0u~ + sin p sin OUy + cos (j)u ~). (19) 
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These jump conditions, together with the original jump condition Eq. , constitutes a set of lowest 
order jump conditions that is enforced in the MIB method. 

In principle, by differentiating these lowest order jump conditions, we can generate even more jump 
conditions. However by doing this, we could create some higher-order derivatives and cross derivatives, 
whose evaluation often involves larger stencils, which is unstable for constructing high-order numerical 
schemes and unfeasible for truly complex interface geometries. Therefore, the spirit of the original MIB 
method is to use only the lowest order interface jump conditions. This spirit is preserved in the proposed 
MIB-FVM. We only use jump conditions in Eqs. (§, @ and 

Thus, at each intersection of the interface and the mesh line, there are four jump conditions described 
by Eqs. and ( p!^ , which involve the function values and only the first order derivatives. 

This property is very important for making the numerical scheme more stable and efficient in handling 
complex interface geometries. Besides, employing only the lowest derivatives will generate a more-banded 
matrix with a smaller conditional number. The enforcement of these jump conditions determines fictitious 
values to be used in the discretization of differential operators in a given PDE. 

The discretization of the differential operator involves six partial derivatives, namely where 

Xj = z. In dealing with complex interface geometries, some of these partial derivatives could be very 
difficult or even impossible to be approximated with some given order of accuracy. Since in constructing a 
second order numerical scheme, there is a pair of fictitious values along each mesh line and they are solved 
simultaneously, with the aforementioned four jump conditions, it is affordable for us to eliminate two of 
them. This redundancy gives two more degrees of freedom for us to design efficient and robust second 
order schemes in handling complex interface geometries. The designed numerical scheme systematically 
eliminates the two partial derivatives that are most difficult to discretize at each intersection of the 
interface and the mesh line by using two redundant jump conditions. In summary, at each intersection of 
the interface and the mesh line, there are four remaining partial derivatives that are need to be discretized. 

The elimination process must comply with the following rules: 

• All fictitious values on all irregular points are determined. At the intersection of the interface and 
each mesh line, two fictitious values are determined at the same time. 


• Along each mesh line, two corresponding derivatives must be kept. 

• Among the remaining four partial derivatives, eliminate two that are most difficult to evaluate 
according to the local geometry. 


From Eq. (15), we have 


Ufj 


= P 


Vjy 


( 20 ) 


Therefore, Eqs. (0,0 and ( p^ can be rewritten as 

/ 


[«c] 


= c 


u 

\ u 


( 21 ) 


where 


C = 



Pilin'^ Pi2l3'^ -Pi2/3~ Pi3/3~'~ 

P21 —P21 P22 —P22 P23 

P31 —P31 P32 —P32 P33 



( 22 ) 
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Here pij is the ijth entry of the transformation matrix C and Ci denotes ith row of C. After eliminating 


the Ith and the mth elements of entry of the array {u '^)^, Eq. (21) becomes 


a[Pu^ + h[ur^\ + c[u(^] — (aCi + bC2 + cC^) 


where 


a = C2lC3rn “ C3iC2n 
b = CsiCim — CiiCsrr 

C = CiiC2m — C2lCirr 


/ ut\ 


y 

u+ 

\u7 / 


(23) 


(24) 


In practice, the two remaining jump conditions Eqs. @ and ( [ 2 ^ are used to determine a pair of fictitious 
values at the intersection point of the interface and each of the mesh line at a time. 


3.2 General MIB scheme 

Consider a geometry illustrated in Eigure|^ The domains are denoted by and . The interface 
intersects the kth ^-mesh line at point {xo^Vo^^o) on fho yz plane. Eor discretizing central derivatives, 
two fictitious values, /ij,/c and are to be determined on irregular grid points. 



j "1 j j+1 j+2 


Eigure 6: Smooth interface at cross section x = Xi. The interface intersects with the kth mesh line in the 
y direction at point {xo^yo^^o)^ where jump conditions need to be matched. At a pair of grid points in 
red color, namely (i, j, k) and (i, j + 1, /c), fictitious values are to be determined. Eour grid points on the 
y mesh line are employed to approximate quantities at {xo^yoi ^ 0)1 and six grid points denoted by black 
squares are employed to approximate u~ \(^Xo.yo,Zo)' 

Here, u~^ u'^ and Uy at the intersection point (xo^yo^zo) are easily expressed by interpolation 
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and the CFD scheme from information in and Q , respectively 

U = (iCq , tCo j-j-l, tCo J-I- 2 ) ' (/i,jf,/c:+ 5 '^ij-|-2,/c) 

Vjy — —1:5+ l—1,/c: + 

Uy = j-, ICi j-i-i, ICi J-I-2) • (/i,jf,/c55'^i,jf+ 2 ,/c) 


T 

(25) 

T 

(26) 

T 

(27) 

T 

5 

(28) 


where denote interpolation weights, which is generated by using the standard Lagrange polynomials 
El. The first subscript s represents either interpolation (when 5 = 0) or first order derivative (when 
s = 1) at point {xq^Vo^ ^ 0)1 while the second subscript represents the node index. 

We only need to compute two out of four remaining first order partial derivatives. For instance, if u~ 
and u~ are easily computed, then by choosing s = 1 and t = 5 in Eqs. (23) and (24), we can eliminate 
14+ and 14+. 

To approximate 14+ or u~ , we need three more function values at the intersection points of the auxiliary 
line y = Ho and the mesh line on the y — z plane. Here we only provide a detailed scheme to approximate 
14“, while other derivatives can be approximated in the same manner. Since these points are not on grid 
nodes, they are interpolated by nearby grid points along the y direction. Therefore, six more auxiliary 
points are involved. As illustrated in Figure^ u~[xo^yo^ Zq) can be approximately as 


Uz {Xo,yo,Zo) = (l^;l,/c,l^^l,/c + l,l^^l,/c+2) • 


Woj 

0 

0 


Woj-^l 

0 

0 


W0,j+2 

0 

0 


0 


0 


^0,i+2 

0 


0 

0 

wlj 


0 

0 

^0,j + l 


0 

0 

'^ 0 ,j +2 

(29) 


where U — (^fij^}^^Uijj^i^]^^Uijj^2,k’)'^i,j,k-\-l’)'^i,j-\-l,k-\-l’)'^i,j-\-2,k-\-l^'^i,j,k-\-2’)'^i,j-\-l,k-\-2’)'^i,j-\-2,k-\-2) • Here 
the superscripts on w denotes different sets of FD weights m- Similarly, 14^ can be approximated along 
the auxiliary line y = yo on the xy plane. Then jump conditions Eqs. and ( [^ combined are used 
to solve the pair of fictitious values and at the same time. Consequently, the two fictitious 

values are expressed as a linear combination of 16 function values at nearby grid points and 4 jump 
conditions. 


3.3 MIB scheme for interface with large curvatures 

The above argument is based on a crucial assumption, which is that on each mesh line, there should 
be enough grid nodes (at least two for employing second order standard ED scheme) near the interface 
inside each subdomain so that the jump conditions can be expressed. However, when the curvature of 
the interface is very large, this assumption cannot be guaranteed for all mesh sizes. 

As is illustrated in Eigure fictitious value fi,j-i,k and cannot be solved along y- 

direction, since there is only one grid point inside the interface in the ^-direction. The concept of 
disassociation is introduced in Ref. [69] in order to disassociate the domain extension from the discretiza¬ 
tion thus broaden the applicability of the MIB method to general interface geometry. In the situation of 
Eigurej^ we cannot solve the fictitious values in the ^-direction, however, there is no problem to solve 
the fictitious values in the ^(-direction. Consequently, the fictitious value at an irregular point, regardless 
of the direction in its calculation, can be used for discretization in any direction involving the grid point 
without loss of accuracy. Eor fictitious values that can be obtained in multiple directions, we may com¬ 
pute their values in the most convenient manner and use them for necessary discretization. In practice, in 
order to make the MIB matrix as symmetric and banded as possible, if the fictitious value can be found 
in the given direction, one should use the fictitious value obtained from the corresponding direction and 
avoid using the disassociation technique. It is also worth of mentioning that disassociation technique 
does not reduce the accuracy of the numerical scheme. Eor instance, if fictitious values obtained from 
the ^-direction has 0{h^) accuracy for some integer m. This accuracy of approximation is independent 
of the direction in which the fictitious values are obtained. Here h is the grid size of the uniform mesh. 
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Figure 7: Disassociation type of irregular grid points at cross section x = Xi. Fictitious value need 
to be deployed at point {i^j — 1,/c) and (i,jf + 1,/c) in order to discretize the PDF at these 

points. However, fi,j-i,k and can not be obtained from the ^-direction by the general MIB 

scheme, since there are not enough points in . Luckily, these fictitious value can be obtained from 
the ^-direction. For the above figure, during the discretization process, fictitious values obtained in the 
2;-direction will be utilized for both the ^-direction and the ^-direction discretizations of the PDE. 


Table 1: Numerical errors and convergence orders for a 6-petal flower interface problem (Case 1). 


Hx X Tly 

L'oo ('^) 

Order 

L2{u) 

Order 

40 X 40 

3.790e-5 


1.910e-5 


80 X 80 

4.867e-6 

2.96 

1.827e-6 

3.38 

160 X 160 

8.545e-7 

2.51 

2.275e-7 

2.99 


4 Numerical studies 

In this section, we examine validity and test the performance of the proposed MIB-FVM for solving the 
3D Poisson equation with discontinuous coefficients. To demonstrate the robustness and accuracy of the 
present MIB-FVM, we carry out many numerical experiments with complex interfaces in both 2D and 3D 
settings. For 2D test cases, we start with a 6-petal flower shape interface which is relatively complex. The 
2D jigsaw puzzle-like interface is also studied. For 3D problems, we consider various interface geometries 
including sphere, ellipsoid, cylinder, 5-petal flower and torus. Finally, we investigate our MIB-FVM’s 
ability to deal with low regularity solutions in two more test cases. 

Case 1. We first consider a 2D interface problem. The 2D Poisson equation is solved in domain 
[—1,1] X [—1,1]. The interface of the problem is a 6-petal flower whose formula in polar coordinate 
system is given as r = ^(1 + | sin(66>)). The designed solution in two different domains are given by 

u'^{x,y) = i + sin(x) sin(^), and u~{x,y) =0. (30) 

The coefficients are given by P~^{x,y) = ^~{x,y) = 1. 

Figure [^depicted the numerical solution and the Loo error on a 40 x 40 x 40 mesh. From the plot of 
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Figure 8: The computed solution on an 80 x 80 mesh (Left chart) and the Loo error (Right chart) for 
Case 1. 


Table 2: Numerical errors and convergence orders for an ellipsoidal interface problem (Case 2). 


nx X Uy 

L^oo ('^) 

Order 

L2 {u) 

Order 

40 X 60 

1.490e-3 


2.176e-4 


80 X 120 

2.255e-4 

2.72 

3.654e-5 

2.57 

160 X 240 

2.448e-5 

3.20 

4.067e-6 

3.16 


the Loo error on the right, it is noticed that the highest error appears at the tips of the petals, where the 
curvature of the interface is relatively higher. Table shows the results of the numerical accuracy tests 
on three successively refined meshes. From the table, it can be seen that MIB-FVM obtains the second 
order accuracy in both the Loo error and L 2 error. 

Case 2. We next consider the 2D jigsaw puzzle-like interface problem. The 2D Poisson equation is 
solved in domain [—1,1] x [0, 3]. The interface of the problem is a jigsaw puzzle-like shape whose formula 
is given as 


I x{0) = 0.6 cos(6>) — 0.3 cos(36>) 

[y{0) = 1.5 + 0.7sin((9) -0.07sin(3<9) + 0.2 sin(7(9) ^ ^ 

The solution in two different subdomains is chosen as 

u+(a;,y) = e*( 2 /^+a;^sin(y)), and u~ {x,y) =+ y'^) (32) 

The discontinuous coefficients are given by /3~^{x^y^z) = 1, and l3~{x^y^z) = 10. 

Figure depicted the computed solution and the Loo error on a 40 x 40 x 40 mesh. The numerical 
errors in terms of Lqo norm are collected in Table which also show the second order convergence of the 
proposed MIB-FVM. 

Case 3. Next we investigate a classical spherical interface problem. The 3D Poisson equation is 
solved in domain [—1,1] X [-1,1] X [—1,1]. To specify the spherical interface, we design the level set 
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Figure 9: The computed solution on a 40 x 40 mesh (Left chart) and the Loo error (Right chart) for Case 

2 . 


Table 3: Numerical errors and convergence orders for a spherical interface problem (Case 3(a)). 


Tlx X fly X riz 

Loo ('^) 

Order 

L2{u) 

Order 

20 X 20 X 20 

1.2771e-3 


5.5733e-4 


40 X 40 X 40 

2.8944e-4 

2.14 

1.2612e-4 

2.14 

80 X 80 X 80 

6.4798e-6 

2.15 

2.8204e-5 

2.16 

160 X 160 X 160 

1.1673e-5 

2.47 

4.6283e-6 

2.61 


function 0(x, z) 


(t>{x,y,z)=rl-{x^+y^+z^), (33) 

where the radius ro = 3/4. The interface is given as T = {(x, y, z) = 0, V(x, y, z) G Two subdomains 
are 0+ = {(x, y, ^)|0 > 0, V(x, y, z) G and = {(x, y^ z)\(j) < 0; V(x, y, z) G f^}. The solution in two 
different subdomains is chosen as 

i4^(x, z) = 0, and u~ {x, y, z) = cos{x) cos{y) cos{z). (34) 

The discontinuous coefficients are given by 

• Case 3(a): 


/3+{x,y,z) =8, and /3 {x,y,z) = l. 


(35) 


• Case 3(b): 

/3+(a;,y,2:) = 1+COS (^ + ^ + ^) , and ^-(a;,2/,2) = 1 +3sin (^ + ^ + ^) . (36) 

Tablelists numerical results for the spherical interface problem. The geometry and error are depicted 
in Figure It is seen that errors are very small at the fine mesh. The second order convergence is 
found. 
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Figure 10: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the error (Right chart) 
for Case 3(a). 



Figure 11: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the error (Right chart) 
for Case 3(b). 


Table 4: Numerical errors and convergence orders for a spherical interface problem (Case 3(b)). 


Ux X riy X Hz 

Lqo ('^) 

Order 

L2 {u) 

Order 

20 X 20 X 20 

8.3466e-5 


3.0069e-5 


40 X 40 X 40 

1.9207e-5 

2.12 

6.8250e-6 

2.14 

80 X 80 X 80 

4.9351e-6 

1.96 

1.5630e-6 

2.13 

160 X 160 X 160 

1.5551e-6 

1.67 

3.7873e-7 

2.05 
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Figure 12: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the error (Right chart) 
for Case 4. 


Table 5: Numerical errors and convergence orders for an ellipsoidal interface problem (Case 4). 


fix x Uy X Hz 

-^00 (^) 

Order 

L2 {u) 

Order 

20 X 20 X 20 

2.1881e-2 


6.8115e-3 


40 X 40 X 40 

5.5719e-3 

1.97 

1.7463e-3 

1.96 

80 X 80 X 80 

1.4675e-3 

1.92 

4.5309e-4 

1.95 

160 X 160 X 160 

3.6068e-4 

2.02 

1.0974e-4 

2.05 


To further investigate the present method, we change the coefficient into a position dependent function 
in Case 3(b). Result is given in Table which also demonstrates the second order accuracy in both L^o 
and 1/2 norms for the solution. 

Case 4. We next consider an ellipsoidal interface problem. The 3D Poisson equation is solved in 
domain [—5,5] x [—5,5] x [—5,5]. The interface of the problem is an ellipsoid given as + (^) ^ 

= 1. To specify the ellipsoidal interface, we define the level set function 0(x,^, z) 

4>{x,y,z) = l-i^(^ ^ (l") ^ (w) ) 

The interface is given as T = {{x^y^z)\(j) = 0,V(x,^,2:) G ft}. Two subdomains are = {{x^y, z)\(l) > 
0,V(x,^, z) G 11}, and = {{x,y,z)\(l) < 0;V(x,^, 2:) G H}. The solution in two different subdomains is 
chosen as 

u^{x^y^z) = e 20 ^ and u~ {x^y, z) = cos{x) cos{y) cos{z). (38) 

The discontinuous coefficients are given by /d+(x, y^z) = z 15, and P~{x, y^ z) = + 10. 

Figure depicted the numerical solution and the Loo error on a 40 x 40 x 40 mesh. From the result 
demonstrated in Table it can be seen that the MIB-FVM obtains the second order accuracy. 

Case 5. In Case 5, the 3D Poisson equation is solved in domain [—4,4] x [—4,4] x [—2,8.4]. The 
interface of the problem is a cylinder of height 27r and base radius tt. The designed solution in two 
different subdomains is chosen as 

{x^y^ z) = X + y ^ z^ and u~{x^y^ z) = co^{x) co^{y) co^{z). (39) 
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Figure 13: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the Loo error (Right chart) 
for Case 5. 


Table 6: Numerical errors and convergence orders for a cylinder interface problem (Case 5). 


Ux X riy X Hz 

Lqo ('^) 

Order 

L2 {u) 

Order 

20 X 20 X 20 

1.2274e-2 


2.9195e-3 


40 X 40 X 40 

3.0505e-3 

2.01 

7.7935e-4 

1.91 

80 X 80 X 80 

7.9267e-4 

1.94 

2.0674e-4 

1.91 

160 X 160 X 160 

1.8598e-4 

2.09 

4.6819e-5 

2.14 
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Figure 14: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the L^o error (Right chart) 
for Case 6. 


Table 7: Numerical errors and convergence orders for a 5-petal flower-base cylinder interface problem 
(Case 6). 


fix X fly X Hz 

Too (^) 

Order 

L2{u) 

Order 

20 X 20 X 20 

4.9246e-2 


1.4819e-3 


40 X 40 X 40 

1.4307e-2 

1.78 

3.1350e-4 

2.24 

80 X 80 X 80 

4.9600e-4 

4.85 

5.6857e-5 

2.46 

160 X 160 X 160 

1.2096e-4 

2.04 

1.5637e-5 

1.86 


The discontinuous coefficients are given by /d+(x, y, z) = 8, and z) = 1. 

Figure p!3| depicted the numerical solution and the Loo error on a 40 x 40 x 40 mesh. From the result 
demonstrated in Table it can be seen that the proposed MIB-FVM achieves the second order accuracy. 

Case 6. In Case 6, the 3D Poisson equation is solved in domain [—5, 5] x [—5, 5] x [—2, 2]. The interface 
of the problem is a 5-petal flower-base cylinder whose polar equation is given by r = | -h | sin(5^) and 
The designed solution in two different subdomains is chosen as 

u~^{x^ y, z) = sm{x) sm{y) sm{z), and u~ {x^y^ z) = . (40) 

The discontinuous coefficients are given by /d^(x, y^ z) = 8, and P~{x^ y^ z) = 1. 

Figure [T^ depicted the numerical solution and the error on a 40 x 40 x 40 mesh. From the result 
demonstrated in Table it can be seen that the present MIB-FVM obtains the second order accuracy. 

Case 7. In Case 7, the 3D Poisson equation is solved in domain [—5,5] x [-5,5] X [—2,2]. The 
interface of the problem is a torus whose equation is given by (3 — x‘^ y‘^Y + = 1- The designed 

solution in two different subdomains is chosen as 

i4^(x,z) = sin(x) sin(^) sin(^) + I, and u~{x^y^ z) = co^{x) co^{y) co^{z). (41) 

The discontinuous coefficients are given by /d+(x, z) = 8, and P~{x^ y^ z) = 1. 

Figure [^depicted the numerical solution and the L^o error on a 40 x 40 x 40 mesh. From the result 
demonstrated in Table it can be seen that the present MIB-FVM attains the second order accuracy. 
Case 8. in this case, we test the proposed method for problems with low solution regularity. 

• Case 8(a): The 3D Poisson equation is solved in domain [—1,1.05] x [—1,1.05] x [—1,1.05]. The 
interface of the problem is a sphere the same as that in Case 3. The designed solution in two 
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Figure 15: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the Loo error (Right chart) 
for Case 7. 


Table 8: Numerical errors and convergence orders for a torus interface problem (Case 7). 


fix X riy X Hz 

-^oo (^) 

Order 

L2 {u) 

Order 

20 X 20 X 20 

3.4040e-2 


4.6519e-3 


40 X 40 X 40 

5.7226e-3 

2.57 

1.0688e-3 

2.12 

80 X 80 X 80 

1.6415e-3 

1.80 

2.6468e-4 

2.01 

160 X 160 X 160 

3.8950e-4 

2.08 

6.9023e-5 

1.94 
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Figure 16: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the L^o error (Right chart) 
for Case 8(a). 



Table 9: Numerical errors and convergence orders for a sphere interface problem (Case 8(a)). 


Tlx X riy X Hz 

Loo ('^) 

Order 

L2 (u) 

Order 

20 X 20 X 20 

1.1578e-3 


3.3263e-4 


40 X 40 X 40 

3.0664e-4 

1.92 

9.2064e-5 

1.85 

80 X 80 X 80 

6.8400e-5 

2.17 

2.0947e-5 

2.14 

160 X 160 X 160 

1.3566e-5 

2.34 

4.4747e-6 

2.24 


different subdomains is chosen as 

u'^{x,y,z) = S, and u~{x,y, z) = {x‘^ + y‘^ ^ z‘^)^. (42) 

The discontinuous coefficients are given by P^{x^y^z) = 4, and P~{x^y^z) = 1. 

• Case 8(b): The 3D Poisson equation is solved in domain [—5,5.05] x [—5,5.05] x [—5,5.05]. The 
interface of the problem is an ellipsoid the same as that in Case 4. The designed solution in two 
different subdomains is chosen as 

i4^(x, z) = 8, and u~{x^ y^ z) = {x‘^ + y‘^-\-z‘^)i^\n.{x ^ yz). (43) 

The discontinuous coefficients are given by /d+(x, y^z) = 2^ + 5, and P~{x, y, z) = + 10. 

In both cases, the second order derivatives of U- blow up at the origin, making the solution only 
continuous in Since the differential form, Eq. Q, is of second order for the diffusion term, 

mathematically, Eq. 0 is invalid even though the underlying conservation law is still valid. However, it 
is noticed that the equivalent term in the integral form, Eq. (|^, involves only the first order derivative. 
This reduction of the derivative order is important in dealing with solution which changes so rapidly in 
space that the spatial derivative does not exist. Eor these reasons, the finite volume method is preferred 
over the finite difference method in solving problems whose solution has less regularity and exhibits local 
discontinuities. 


21 













Figure 17: The computed solution on a 40 x 40 x 40 mesh (Left chart) and the Loo error (Right chart) 
for Case 8(b). 


Table 10: Numerical errors and convergence orders for an ellipsoid interface problem (Case 8(b)). 


Tlx X riy X riz 

AqO ('^) 

Order 

L2{u) 

Order 

20 X 20 X 20 

3.2934e-2 


I.I797e-2 


40 X 40 X 40 

7.5420e-3 

2.13 

2.9620e-3 

1.99 

80 X 80 X 80 

I.9I83e-3 

1.98 

7.5382e-4 

1.97 

160 X 160 X 160 

4.7448e-4 

2.02 

I.870Ie-4 

2.01 


In Case 8(a), we choose constant diffusive coefficients /3+ and whereas in Case 8(b), and [3~ 
are both position dependent. The numerical solutions and the Loo errors of Case 8(a) and Case 8(b) are 
depicted in Figure and Figure pT| respectively. From the numerical results collected in Table and 
Table pLj it is seen that the second order accuracy in both cases are essentially obtained. 


5 Conclusion 

In this work, we introduce the finite volume formulation of matched interface and boundary method 
(MIB-FVM) for solving elliptic interface problems arising from modeling material interface in practice. 
Much effort has been taken to develop advanced numerical schemes for such problems in the past few 
decades. However, challenges still remain in this field. One of the challenge concerns the development of 
methods with higher order accuracy. Another challenge is to develop methods for dealing with complex 
interface geometries. The matched interface and boundary (MIB) method has been proved to be able 
to deal with these challenges. However, since the MIB method is based on the collocation formulation, 
it does not work well for solutions with low regularity. Based on the integral form rather than the 
differential form, the finite volume method (FVM) is known for being able to deal with problems with 
low regularity solutions and better conserve mass and flux. Second order finite volume methods have also 
been formulated in literature for elliptic interface problems. Motivated by these successes, we propose 
the present MIB-FVM to take the advantages of MIB and FVM. 

To reduce the computational cost of mesh generation, we utilize the Cartesian mesh, although the 
proposed MIB-FVM can be realized on irregular meshes as well. We also employ the vertex-centered 
FVM, for which the cubic control volumes are associated with the grid point at the center. The control 
volumes are categorized into regular and irregular types, where special treatments of the MIB formulation 
are needed for the irregular ones in order to maintain the designed accuracy. We study the proposed MIB- 
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FVM by a number of classical test cases, including both 2D cases such as 6-petal flower and jigsaw-puzzle 
like shape, and 3D cases such as sphere, ellipsoid, standard cylinder, flower-based cylinder and torus. 
The numerical results all demonstrate the second order accuracy of the proposed method. The ability of 
handling interface with complex geometries and solutions with low regularity continuous) indicates 
that the proposed MIB-FVM combines the advantage of both MIB and FVM in meeting numerical 
challenges. Future work will be done to improve the presnt method in dealing with even more complex 
interface geometries, for example, interface with geometric singularities which are commonly seen in 
practical applications such as electrostatic analysis of proteins. 
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